Mixture model for blocked ring canals
Jonathan Harrison
24/09/2018
Blocking of Ring canals
Forward simulating
Forward simulating from the model suggests that blocking of ring canals could account for differences between observations and observed counts of mRNA particles in cells far from the oocyte.
##
## Gradient evaluation took 1.7e-05 seconds
## 1000 transitions using 10 leapfrog steps per transition would take 0.17 seconds.
## Adjust your expectations accordingly!
##
##
##
## Elapsed Time: 1.58321 seconds (Warm-up)
## 1.57348 seconds (Sampling)
## 3.15668 seconds (Total)
##
##
## Gradient evaluation took 1.7e-05 seconds
## 1000 transitions using 10 leapfrog steps per transition would take 0.17 seconds.
## Adjust your expectations accordingly!
##
##
##
## Elapsed Time: 1.54349 seconds (Warm-up)
## 1.62903 seconds (Sampling)
## 3.17252 seconds (Total)
##
##
## Gradient evaluation took 8e-06 seconds
## 1000 transitions using 10 leapfrog steps per transition would take 0.08 seconds.
## Adjust your expectations accordingly!
##
##
##
## Elapsed Time: 1.52726 seconds (Warm-up)
## 1.59747 seconds (Sampling)
## 3.12473 seconds (Total)
##
##
## Gradient evaluation took 9e-06 seconds
## 1000 transitions using 10 leapfrog steps per transition would take 0.09 seconds.
## Adjust your expectations accordingly!
##
##
##
## Elapsed Time: 1.61221 seconds (Warm-up)
## 1.59713 seconds (Sampling)
## 3.20933 seconds (Total)
Mixture model
Blocking each of the 15 ring canals, or not blocking any, results in different distributions of mRNA across the cells in the tissue at steady state. We consider using this blocking mechanism as a more mechanistic noise model.
We fit a mixture model with 16 hidden states corresponding to which (if any) ring canal is blocked. We marginalise over the hidden states by summing over these states and infer the ring canal transport bias, \(\nu\). The model is as follows: \[\begin{align} p(\nu | \mathbf{y}) &\propto p(\mathbf{y} | \nu) p(\nu) \\ &= \sum_{z=1}^{16} \theta_z p(\mathbf{y} | \nu, z) p(\nu) \\ &= \sum_{z=1}^{16} \theta_z N(\mathbf{y} \, | \, \mathbf{k}_z(\nu)) p(\nu) \\ \end{align}\]where \(\mathbf{y}\) is the observed mRNA counts across cells, \(z \in \{ 1,2, \dots, 16\}\), \(\nu\) is the transport bias, \(\mathbf{k}\) is the steady state distribution of mRNA across cells predicted from the model when blocking the \(z\)th ring canal, and \(\theta_z\) are the mixture weights.
Observed distribution of RNA across cells
To illustrate these ideas, we plot the RNA distributions for WT and OE based on observed data. 
Simulated distribution of RNA across cells
We can generate similar distributions from the mixture model. In the plot below, we show the distribution characteristic of each of the 16 states in the mixture corresponding to blocking a different ring canal (or no blocking).
Analysis of ring canal bias at steady state
By considering the distribution of RNA across cells at steady state, we were able to infer the transport bias \(\nu\) for transport through ring canals.
Wild type
## [1] "using real data \n"
## Inference for Stan model: model_comparison_at_stst2.
## 4 chains, each with iter=2000; warmup=1000; thin=1;
## post-warmup draws per chain=1000, total post-warmup draws=4000.
##
## mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
## nu 0.97 0 0.02 0.94 0.96 0.98 0.99 1.00 1892 1
## xi 0.10 0 0.01 0.09 0.10 0.10 0.11 0.11 2288 1
## phi 0.35 0 0.05 0.26 0.32 0.35 0.38 0.45 2460 1
##
## Samples were drawn using NUTS(diag_e) at Thu Oct 4 12:14:21 2018.
## For each parameter, n_eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor on split chains (at
## convergence, Rhat=1).
Posterior predictive distribution: 
Overexpressor
## [1] "using real data \n"
## [1] "swapped"
## [1] 0.00001 2.40820 12.43159 14.79800 19.54926 21.98536 22.69240 22.81189
## [9] 24.29904
## [1] 9
## Inference for Stan model: model_comparison_at_stst2.
## 4 chains, each with iter=2000; warmup=1000; thin=1;
## post-warmup draws per chain=1000, total post-warmup draws=4000.
##
## mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
## nu 0.99 0 0.01 0.95 0.98 0.99 1.00 1.00 3354 1
## xi 0.40 0 0.02 0.36 0.38 0.40 0.41 0.44 3422 1
## phi 0.35 0 0.05 0.26 0.32 0.35 0.38 0.45 2974 1
##
## Samples were drawn using NUTS(diag_e) at Thu Oct 4 12:18:00 2018.
## For each parameter, n_eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor on split chains (at
## convergence, Rhat=1).
Posterior predictive distribution: 
As expected, we obtain a value of \(\nu\) close to 1, indicating strong bias in the direction of the oocyte. The value/distribution of \(\nu\) is similar for both the WT and OE. The main difference between fitting to the WT and OE data is the estimated value for the noise in the OE data. This motivated the use of the different noise model to examine differences between the WT and OE.
Using the mixture model as a measurement error model
In order to simulate data similar to the overexpression data, this model must be dominated by a large value of the noise. The measurement noise we use is not very mechanistic - simply Gaussian noise here. Since we have proposed a biological mechanism leading to this noise, we instead consider using this to motivate a different measurement error model. We use a mixture model and sum over all the possible blocked ring canals to enable to infer the transport bias \(\nu\) as usual.
Overexpressor with mixture noise model for blocked ring canals
## [1] "using real data \n"
## [1] "swapped"
## [1] 0.00001 2.40820 12.43159 14.79800 19.54926 21.98536 22.69240 22.81189
## [9] 24.29904
## [1] 9
## Inference for Stan model: model_comparison_at_stst3.
## 4 chains, each with iter=2000; warmup=1000; thin=1;
## post-warmup draws per chain=1000, total post-warmup draws=4000.
##
## mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
## nu 0.48 0 0.06 0.39 0.45 0.48 0.52 0.61 2751 1
## xi 0.31 0 0.03 0.25 0.29 0.31 0.33 0.38 4000 1
## phi 0.34 0 0.03 0.27 0.32 0.35 0.37 0.41 4000 1
##
## Samples were drawn using NUTS(diag_e) at Thu Oct 4 13:11:58 2018.
## For each parameter, n_eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor on split chains (at
## convergence, Rhat=1).
Model comparison to compare blocking to no blocking
The mixture model at steady state seems to fit reasonably well. For some cells such as cell 2, the credible interval is much wider than the observed data. Is the mixture model (describing the blocking mechanism) genuinely an improvement on the previous model without blocking at steady state? We use model comparison via leave one out cross validation to find out.
##
## Computed from 4000 by 9 log-likelihood matrix
##
## Estimate SE
## elpd_loo -45.2 14.8
## p_loo 2.3 1.0
## looic 90.3 29.5
## ------
## Monte Carlo SE of elpd_loo is 0.0.
##
## All Pareto k estimates are good (k < 0.5).
## See help('pareto-k-diagnostic') for details.
##
## Computed from 4000 by 9 log-likelihood matrix
##
## Estimate SE
## elpd_loo -35.9 7.9
## p_loo 4.3 1.4
## looic 71.7 15.8
## ------
## Monte Carlo SE of elpd_loo is 0.1.
##
## Pareto k diagnostic values:
## Count Pct. Min. n_eff
## (-Inf, 0.5] (good) 8 88.9% 1243
## (0.5, 0.7] (ok) 1 11.1% 262
## (0.7, 1] (bad) 0 0.0% <NA>
## (1, Inf) (very bad) 0 0.0% <NA>
##
## All Pareto k estimates are ok (k < 0.7).
## See help('pareto-k-diagnostic') for details.## elpd_diff se
## 9.3 8.0Some extra models: